Escaping Earth with Julia

Copyright 2016, Pedro Belin Castellucci,

This work is licensed under a Creative Commons Attribution 4.0 International License.</i>

Gravity is a familiar concept. Every object near the Earth is subject to Earth gravitational force, which attract the object to its center. This gravitational force is associated with gravitational acceleration $g$. Near the surface of the Earth $g \approx 9.8\ m/s²$. However, assuming the acceleration is constant may be unrealistic for some applications. In fact, different parameters affect $g$ such as latitude and altitude. In this Notebook, we will play a little bit with the dependence of the gravitational acceleration and altitude. Our goal is not to do rigorous Math or Physics, but to use the context to learn some Julia programming basics (for some Math and Physics, the reader is refered to this link).

Instead of assuming a constant acceleration, a more realistic model for the gravity acceleration $g_e$ is:

$$g_e = g_0 \Big(\frac{r_e}{r_e + h}\Big)^2,$$

in which $g_0$ is called the standard acceleration of gravity and is defined as $g_0 = 9.80665\ m/s^2$. The actual acceleration is $g_e$ and depends on $r_e$ and $h$, the average radius of Earth and the height (above sea level) of the point from which we are computing $g_e$, respectively.

Defining variables

Let us define Julia variables to keep $g_0 = 9.80665\ m/s^2$ and the average radius of earth $r_e = 6 378 000\ m$.


In [1]:
g0 = 9.80665  # g (m/s²)
re = 6378000  # Average radius of Earth in meters (m).


Out[1]:
6378000

Note that we do not need to declare the type of our variables. However they do have a type, which we can check with typeof.


In [2]:
typeof(g0), typeof(re), typeof("Earth")


Out[2]:
(Float64,Int64,String)

Julia supports Unicode variable names. This allows the definition of variables with more mathy names, e.g. with subscripts or with greek symbols. To declare a variable $\alpha$, for example, we can type \alpha{TAB}. Thus, let us redefine our variables with the g\_0{TAB} and r\_e{TAB}.


In [3]:
g₀ = g0
rₑ = re


Out[3]:
6378000

Defining a function

Next, we define a function called localGravity to compute $g_e$ as function of $g_0$, $r_e$ and $h$.


In [4]:
function localGravity(g₀, rₑ, h)
    g₀*(rₑ/(rₑ + h))^2  # This is the return of the function
end


Out[4]:
localGravity (generic function with 1 method)

Note that the return of a function is the statement in its last line. Another way of defining the same function is using a one line definition.


In [5]:
localGravity(g₀, rₑ, h) = g₀*(rₑ/(rₑ + h))^2


WARNING: Method definition localGravity(Any, Any, Any) in module Main at In[4]:2 overwritten at In[5]:1.
Out[5]:
localGravity (generic function with 1 method)

We get a warning because we are redefining the localGravity function.

While-loop

We will use a while-loop to simulate the position of our particle over time. As a first example, let us just see how to print out the time instants from 0 to 4s with a time step of $\Delta t=0.5s$.


In [6]:
Δt = 0.5  # Time step.
t = 0  # Initial time.
T = 4  # Final time.

while t <= T
    println(t, "s")
    t += Δt
end


0s
0.5s
1.0s
1.5s
2.0s
2.5s
3.0s
3.5s
4.0s

Using a while-loop, we can simulate a vertical launch of a particle. We will assume that the variation of position and velocity of our particle is given by:

$$h_k = h_{k-1} + v_{k-1} \Delta t,$$$$v_k = v_{k-1} - g_e^{k-1} \Delta t,$$

in which $h_k$ and $v_k$ are the height and velocity at time $k$, respectively. $g_e^{k-1}$ is given by the localGravity function and $\Delta t$ is the discrete time step we choose for our simulation. As a first example, we will consider $(h_0, v_0) = (0, 10)$ and simulate 2.2 seconds with $\Delta t = 0.1s.$


In [7]:
# Initial conditions:
v0 = 10.0  # m/s
h0 = 0.0   # m
time = 0.0  # s

# Simulation time:
endOfSimulation = 2.2

# Time step:
Δt = 0.1  # s

while time <= endOfSimulation

    # Printing a formatted output:
    toPrint = @sprintf("t = %.2fs \t h = %.2fm", time, h0)
    
    println(toPrint)
    # Updating height and velocity:
    
    h0, v0 = h0 + v0*Δt, v0 - localGravity(g₀, rₑ, h0)*Δt
    
    # Incrementing time:
    time += Δt
end


t = 0.00s 	 h = 0.00m
t = 0.10s 	 h = 1.00m
t = 0.20s 	 h = 1.90m
t = 0.30s 	 h = 2.71m
t = 0.40s 	 h = 3.41m
t = 0.50s 	 h = 4.02m
t = 0.60s 	 h = 4.53m
t = 0.70s 	 h = 4.94m
t = 0.80s 	 h = 5.25m
t = 0.90s 	 h = 5.47m
t = 1.00s 	 h = 5.59m
t = 1.10s 	 h = 5.61m
t = 1.20s 	 h = 5.53m
t = 1.30s 	 h = 5.35m
t = 1.40s 	 h = 5.08m
t = 1.50s 	 h = 4.70m
t = 1.60s 	 h = 4.23m
t = 1.70s 	 h = 3.66m
t = 1.80s 	 h = 3.00m
t = 1.90s 	 h = 2.23m
t = 2.00s 	 h = 1.37m
t = 2.10s 	 h = 0.41m

We printed out the height of the particle as time evolves. To make it a little prettier, we used the macro @sprintf. We are not going to discuss macros here, but superficially they are like functions.

However, we can do a better job with the output. Julia has different plotting libraries, we are going to use Plots to plot the trajectory of the particle over time. To install it you can use the following:


In [8]:
Pkg.add("Plots");


INFO: Cloning cache of ColorTypes from https://github.com/JuliaGraphics/ColorTypes.jl.git
INFO: Cloning cache of Colors from https://github.com/JuliaGraphics/Colors.jl.git
INFO: Cloning cache of DataStructures from https://github.com/JuliaCollections/DataStructures.jl.git
INFO: Cloning cache of FixedPointNumbers from https://github.com/JuliaMath/FixedPointNumbers.jl.git
INFO: Cloning cache of FixedSizeArrays from https://github.com/SimonDanisch/FixedSizeArrays.jl.git
INFO: Cloning cache of Measures from https://github.com/JuliaGraphics/Measures.jl.git
INFO: Cloning cache of NaNMath from https://github.com/mlubin/NaNMath.jl.git
INFO: Cloning cache of PlotThemes from https://github.com/JuliaPlots/PlotThemes.jl.git
INFO: Cloning cache of PlotUtils from https://github.com/JuliaPlots/PlotUtils.jl.git
INFO: Cloning cache of Plots from https://github.com/JuliaPlots/Plots.jl.git
INFO: Cloning cache of RecipesBase from https://github.com/JuliaPlots/RecipesBase.jl.git
INFO: Cloning cache of Reexport from https://github.com/simonster/Reexport.jl.git
INFO: Cloning cache of Showoff from https://github.com/JuliaGraphics/Showoff.jl.git
INFO: Cloning cache of SpecialFunctions from https://github.com/JuliaMath/SpecialFunctions.jl.git
INFO: Cloning cache of StatsBase from https://github.com/JuliaStats/StatsBase.jl.git
INFO: Installing ColorTypes v0.5.2
INFO: Installing Colors v0.7.4
INFO: Installing DataStructures v0.6.1
INFO: Installing FixedPointNumbers v0.3.9
INFO: Installing FixedSizeArrays v0.2.5
INFO: Installing Measures v0.1.0
INFO: Installing NaNMath v0.2.6
INFO: Installing PlotThemes v0.1.4
INFO: Installing PlotUtils v0.4.4
INFO: Installing Plots v0.11.4
INFO: Installing RecipesBase v0.1.0
INFO: Installing Reexport v0.1.0
INFO: Installing Showoff v0.1.1
INFO: Installing SpecialFunctions v0.2.0
INFO: Installing StatsBase v0.17.0
INFO: Building Plots
INFO: Cannot find deps/plotly-latest.min.js... downloading latest version.
--2018-01-27 16:48:01--  https://cdn.plot.ly/plotly-latest.min.js
Resolving cdn.plot.ly (cdn.plot.ly)... 151.101.52.69
Connecting to cdn.plot.ly (cdn.plot.ly)|151.101.52.69|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 2413253 (2.3M) [application/javascript]
Saving to: '/home/pedro/.julia/v0.5/Plots/deps/plotly-latest.min.js’

     0K .......... .......... .......... .......... ..........  2%  112K 21s
    50K .......... .......... .......... .......... ..........  4%  307K 14s
   100K .......... .......... .......... .......... ..........  6%  327K 11s
   150K .......... .......... .......... .......... ..........  8%  342K 10s
   200K .......... .......... .......... .......... .......... 10%  349K 9s
   250K .......... .......... .......... .......... .......... 12%  375K 8s
   300K .......... .......... .......... .......... .......... 14%  357K 8s
   350K .......... .......... .......... .......... .......... 16%  290K 7s
   400K .......... .......... .......... .......... .......... 19%  399K 7s
   450K .......... .......... .......... .......... .......... 21%  450K 6s
   500K .......... .......... .......... .......... .......... 23%  409K 6s
   550K .......... .......... .......... .......... .......... 25%  454K 6s
   600K .......... .......... .......... .......... .......... 27%  447K 5s
   650K .......... .......... .......... .......... .......... 29%  482K 5s
   700K .......... .......... .......... .......... .......... 31%  483K 5s
   750K .......... .......... .......... .......... .......... 33%  367K 5s
   800K .......... .......... .......... .......... .......... 36%  513K 4s
   850K .......... .......... .......... .......... .......... 38%  521K 4s
   900K .......... .......... .......... .......... .......... 40%  512K 4s
   950K .......... .......... .......... .......... .......... 42%  612K 4s
  1000K .......... .......... .......... .......... .......... 44%  494K 4s
  1050K .......... .......... .......... .......... .......... 46%  650K 3s
  1100K .......... .......... .......... .......... .......... 48%  525K 3s
  1150K .......... .......... .......... .......... .......... 50%  500K 3s
  1200K .......... .......... .......... .......... .......... 53%  578K 3s
  1250K .......... .......... .......... .......... .......... 55%  685K 3s
  1300K .......... .......... .......... .......... .......... 57%  660K 3s
  1350K .......... .......... .......... .......... .......... 59%  630K 2s
  1400K .......... .......... .......... .......... .......... 61%  703K 2s
  1450K .......... .......... .......... .......... .......... 63%  607K 2s
  1500K .......... .......... .......... .......... .......... 65%  768K 2s
  1550K .......... .......... .......... .......... .......... 67%  492K 2s
  1600K .......... .......... .......... .......... .......... 70%  805K 2s
  1650K .......... .......... .......... .......... .......... 72%  766K 2s
  1700K .......... .......... .......... .......... .......... 74%  688K 1s
  1750K .......... .......... .......... .......... .......... 76%  839K 1s
  1800K .......... .......... .......... .......... .......... 78%  765K 1s
  1850K .......... .......... .......... .......... .......... 80%  768K 1s
  1900K .......... .......... .......... .......... .......... 82%  846K 1s
  1950K .......... .......... .......... .......... .......... 84%  609K 1s
  2000K .......... .......... .......... .......... .......... 86%  788K 1s
  2050K .......... .......... .......... .......... .......... 89%  952K 1s
  2100K .......... .......... .......... .......... .......... 91%  877K 0s
  2150K .......... .......... .......... .......... .......... 93%  805K 0s
  2200K .......... .......... .......... .......... .......... 95%  932K 0s
  2250K .......... .......... .......... .......... .......... 97%  965K 0s
  2300K .......... .......... .......... .......... .......... 99%  844K 0s
  2350K ......                                                100%  826K=4.7s

2018-01-27 16:48:07 (497 KB/s) - '/home/pedro/.julia/v0.5/Plots/deps/plotly-latest.min.js’ saved [2413253/2413253]

INFO: Package database updated
INFO: METADATA is out-of-date — you may not have the latest version of Plots
INFO: Use `Pkg.update()` to get the latest versions of your packages

Ranges, for-loops, conditions and arrays

Before using the Plots library, we will take a look at ranges, for-loops and conditions. Ranges are objects that act like arrays in most cases but are implemented in an efficient way. Let us create a range.


In [9]:
myRange = 1:10


Out[9]:
1:10

We can check the type of the object.


In [10]:
typeof(myRange)


Out[10]:
UnitRange{Int64}

To check the elements inside the range object, we can use collect.


In [11]:
collect(myRange)


Out[11]:
10-element Array{Int64,1}:
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10

Note that the syntax a:b creates a collection of the integer number in the interval $[a, b]$. By default, each element is separated by a unit but we can also specify the step of the range to be created.


In [12]:
collect(1:0.1:2)


Out[12]:
11-element Array{Float64,1}:
 1.0
 1.1
 1.2
 1.3
 1.4
 1.5
 1.6
 1.7
 1.8
 1.9
 2.0

Exercise: Use ranges to create the following arrays:

  1. Even numbers in the interval [10, 40];
  2. Odd numbers in the interval [5, 15];
  3. [-10.5, -11, -11.5, -12];
  4. An empty array;

In [ ]:

For-loops are suitable to iterate over range objects.


In [14]:
for i in myRange
    print("$i ") 
end


1 2 3 4 5 6 7 8 9 10 

In Julia, true and false are Bool types.


In [15]:
5 == 2, 5 >= 2, 1 < 2 && 2 < 3, 1 < 2 < 3, typeof(5 == 2)


Out[15]:
(false,true,true,true,Bool)

The Bool types can be used in conditional statements like the one that follows.


In [16]:
a, b = 5, 2
if a == b
    println("The numbers are equal.")
else
    println("The numbers are different.")
end


The numbers are different.

Exercise: Write an algorithm to tell which integer numbers between 1 and 10 are even and odd.


In [ ]:

Arrays

In Julia, arrays have two usages. They can be used as lists and as a mathematical vector. Let us declare an array to check these behaviors.


In [17]:
a = [1, 3, 4, 2]


Out[17]:
4-element Array{Int64,1}:
 1
 3
 4
 2

We have an array with four Int64 elements. Arrays are indexed starting at one, knowing that, try to predict the outcome of each following statement.


In [18]:
println(a[1])
println(a[2])
println(a[end])
println(a[1:3])
println(a[1:end])


1
3
2
[1,3,4]
[1,3,4,2]

Exercise: try to create an array with values of different types.


In [19]:
b = [1, 2.0, 3, "String"]


Out[19]:
4-element Array{Any,1}:
 1        
 2.0      
 3        
  "String"

We can add elements to an array using the push! method.


In [20]:
push!(a, 10)


Out[20]:
5-element Array{Int64,1}:
  1
  3
  4
  2
 10

As mentioned, arrays can also be used in a mathematical context.


In [21]:
println(2a)
println(2 + a)
println([1, 2, 3] + [1, 2, 3])


[2,6,8,4,20]
[3,5,6,4,12]
[2,4,6]

Julia also provides different mathematical operations. Like the dot product.


In [22]:
dot([1, 2, 3], [1, 2, 3])


Out[22]:
14

And the cross product.


In [23]:
cross([1, 2, 3], [1, 1, 1])


Out[23]:
3-element Array{Int64,1}:
 -1
  2
 -1

We can check all the operations available issuing the following command.


In [24]:
methodswith(Array)


Out[24]:
237-element Array{Method,1}:
  • +(A::Array, B::SparseMatrixCSC) at sparse/sparsematrix.jl:1711
  • +(A::SparseMatrixCSC, B::Array) at sparse/sparsematrix.jl:1709
  • -(A::SparseMatrixCSC, B::Array) at sparse/sparsematrix.jl:1714
  • -(A::Array, B::SparseMatrixCSC) at sparse/sparsematrix.jl:1716
  • .*{T<:AbstractString}(v::Array{T,1}, s::AbstractString) at strings/basic.jl:85
  • .*{T<:AbstractString}(s::AbstractString, v::Array{T,1}) at strings/basic.jl:86
  • ./(A::SparseMatrixCSC, B::Array) at sparse/sparsematrix.jl:1724
  • ./(A::Array, B::SparseMatrixCSC) at sparse/sparsematrix.jl:1725
  • .\(A::SparseMatrixCSC, B::Array) at sparse/sparsematrix.jl:1730
  • .\(A::Array, B::SparseMatrixCSC) at sparse/sparsematrix.jl:1731
  • .^{T<:Integer}(A::BitArray, B::Array{T,N<:Any}) at broadcast.jl:406
  • .^(A::SparseMatrixCSC, B::Array) at sparse/sparsematrix.jl:1739
  • .^(A::Array, B::SparseMatrixCSC) at sparse/sparsematrix.jl:1740
  • A_ldiv_B!{Tb<:Complex{T<:Real}}(x::Array{Tb,1}, lu::Base.SparseArrays.UMFPACK.UmfpackLU{Float64,Ti<:Union{Int32,Int64}}, b::Array{Tb,1}) at sparse/umfpack.jl:394
  • A_ldiv_B!{T<:Union{Complex{Float64},Float64}}(lu::Base.SparseArrays.UMFPACK.UmfpackLU{T,Ti<:Union{Int32,Int64}}, b::Array{T,1}) at sparse/umfpack.jl:390
  • A_ldiv_B!{T<:Union{Complex{Float64},Float64}}(lu::Base.SparseArrays.UMFPACK.UmfpackLU{T,Ti<:Union{Int32,Int64}}, b::Array{T,2}) at sparse/umfpack.jl:391
  • A_ldiv_B!{Tb<:Complex{T<:Real}}(lu::Base.SparseArrays.UMFPACK.UmfpackLU{Float64,Ti<:Union{Int32,Int64}}, b::Array{Tb,1}) at sparse/umfpack.jl:406
  • Ac_ldiv_B!{Tb<:Complex{T<:Real}}(x::Array{Tb,1}, lu::Base.SparseArrays.UMFPACK.UmfpackLU{Float64,Ti<:Union{Int32,Int64}}, b::Array{Tb,1}) at sparse/umfpack.jl:394
  • Ac_ldiv_B!{T<:Union{Complex{Float64},Float64}}(lu::Base.SparseArrays.UMFPACK.UmfpackLU{T,Ti<:Union{Int32,Int64}}, b::Array{T,1}) at sparse/umfpack.jl:390
  • Ac_ldiv_B!{T<:Union{Complex{Float64},Float64}}(lu::Base.SparseArrays.UMFPACK.UmfpackLU{T,Ti<:Union{Int32,Int64}}, b::Array{T,2}) at sparse/umfpack.jl:391
  • Ac_ldiv_B!{Tb<:Complex{T<:Real}}(lu::Base.SparseArrays.UMFPACK.UmfpackLU{Float64,Ti<:Union{Int32,Int64}}, b::Array{Tb,1}) at sparse/umfpack.jl:406
  • At_ldiv_B!{Tb<:Complex{T<:Real}}(x::Array{Tb,1}, lu::Base.SparseArrays.UMFPACK.UmfpackLU{Float64,Ti<:Union{Int32,Int64}}, b::Array{Tb,1}) at sparse/umfpack.jl:394
  • At_ldiv_B!{T<:Union{Complex{Float64},Float64}}(lu::Base.SparseArrays.UMFPACK.UmfpackLU{T,Ti<:Union{Int32,Int64}}, b::Array{T,1}) at sparse/umfpack.jl:390
  • At_ldiv_B!{T<:Union{Complex{Float64},Float64}}(lu::Base.SparseArrays.UMFPACK.UmfpackLU{T,Ti<:Union{Int32,Int64}}, b::Array{T,2}) at sparse/umfpack.jl:391
  • At_ldiv_B!{Tb<:Complex{T<:Real}}(lu::Base.SparseArrays.UMFPACK.UmfpackLU{Float64,Ti<:Union{Int32,Int64}}, b::Array{Tb,1}) at sparse/umfpack.jl:406
  • PipeBuffer(data::Array{UInt8,1}, maxsize::Int64) at iobuffer.jl:33
  • PipeBuffer(data::Array{UInt8,1}) at iobuffer.jl:33
  • abs(f::Array{Base.Pkg.Resolve.MaxSum.FieldValues.FieldValue,1}) at pkg/resolve/fieldvalue.jl:67
  • append!(A::Array{Bool,1}, items::BitArray{1}) at bitarray.jl:671
  • append!{T}(a::Array{T,1}, items::AbstractArray{T<:Any,1}) at array.jl:492
  • ascii(v::Array{UInt8,1}) at deprecated.jl:49
  • bytestring(v::Array{UInt8,1}) at deprecated.jl:49
  • convert(::Type{String}, v::Array{UInt8,1}) at strings/string.jl:250
  • convert(::Type{String}, a::Array{UInt8,1}, invalids_as::AbstractString) at deprecated.jl:49
  • convert{T,n}(::Type{Array{T,n}}, x::Array{T,n}) at array.jl:233
  • convert{T,n}(::Type{Array{T,N<:Any}}, x::Array{T,n}) at array.jl:232
  • convert{N}(::Type{BitArray{N}}, A::Array{Bool,N}) at bitarray.jl:527
  • convert{TS,TA,N}(::Type{SharedArray{TS,N}}, A::Array{TA,N}) at sharedarray.jl:280
  • convert{T}(::Type{SharedArray{T,N<:Any}}, A::Array) at sharedarray.jl:279
  • convert(::Type{SharedArray}, A::Array) at sharedarray.jl:278
  • copy{T<:Array{T,N}}(a::T) at array.jl:70
  • copy!{T}(dest::Array{T,N<:Any}, doffs::Integer, src::Array{T,N<:Any}, soffs::Integer, n::Integer) at array.jl:60
  • copy!(dest::BitArray, doffs::Integer, src::Array, soffs::Integer, n::Integer) at bitarray.jl:462
  • copy!{T}(dest::Array{T,N<:Any}, src::Array{T,N<:Any}) at array.jl:68
  • copy!(dest::BitArray, src::Array) at bitarray.jl:471
  • copy!(S::SharedArray, A::Array) at sharedarray.jl:454
  • copy!{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},Ti<:Integer}(dest::Array{T,N<:Any}, rdest::Union{Range{Ti},UnitRange{Ti}}, src::Array{T,N<:Any}, rsrc::Union{Range{Ti},UnitRange{Ti}}) at linalg/blas.jl:1403
  • display(d::Base.REPL.REPLDisplay, md::Array{Base.Markdown.MD,1}) at markdown/Markdown.jl:62
  • div(A::BitArray, B::Array{Bool,N<:Any}) at bitarray.jl:1112
  • div(A::Array{Bool,N<:Any}, B::BitArray) at bitarray.jl:1113
  • done(a::Array, i) at array.jl:381
  • dot{T<:Union{Float32,Float64},TI<:Integer}(x::Array{T,1}, rx::Union{Range{TI},UnitRange{TI}}, y::Array{T,1}, ry::Union{Range{TI},UnitRange{TI}}) at linalg/matmul.jl:48
  • dot{T<:Union{Complex{Float32},Complex{Float64}},TI<:Integer}(x::Array{T,1}, rx::Union{Range{TI},UnitRange{TI}}, y::Array{T,1}, ry::Union{Range{TI},UnitRange{TI}}) at linalg/matmul.jl:61
  • dump(io::IO, x::Array, n::Int64, indent) at show.jl:1132
  • eigvecs{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},Eigenvalue<:Real}(A::SymTridiagonal{T}, eigvals::Array{Eigenvalue,1}) at linalg/tridiag.jl:179
  • evalfile(path::AbstractString, args::Array{String,1}) at loading.jl:504
  • fill!{T<:Union{AbstractFloat,Integer}}(a::Array{T,N<:Any}, x) at array.jl:157
  • flipdim{T}(A::Array{T,N<:Any}, d::Integer) at arraymath.jl:137
  • frexp{T<:AbstractFloat}(A::Array{T,N<:Any}) at math.jl:333
  • gensym(a::Array{UInt8,1}) at expr.jl:8
  • getindex(A::Array, i1::Real) at array.jl:386
  • getindex(A::Array, i1::Real, i2::Real, I::Real...) at array.jl:387
  • getindex(A::Array, I::UnitRange{Int64}) at array.jl:391
  • getindex(A::Array, c::Colon) at array.jl:401
  • getindex{S,T<:Real}(A::Array{S,N<:Any}, I::Range{T}) at array.jl:411
  • indmax(f::Array{Base.Pkg.Resolve.MaxSum.FieldValues.FieldValue,1}) at pkg/resolve/fieldvalue.jl:85
  • insert!{T}(a::Array{T,1}, i::Integer, item) at array.jl:554
  • isassigned{T}(a::Array{T,N<:Any}, i::Int64...) at array.jl:34
  • kill(ps::Array{Base.Process,1}) at process.jl:646
  • launch(manager::Base.SSHManager, params::Dict, launched::Array, launch_ntfy::Condition) at managers.jl:119
  • launch(manager::Base.LocalManager, params::Dict, launched::Array, c::Condition) at managers.jl:312
  • length(a::Array) at array.jl:29
  • lexcmp(a::Array{UInt8,1}, b::Array{UInt8,1}) at array.jl:688
  • map(f::Function, a::Array{Any,1}) at essentials.jl:124
  • merge!(repo::Base.LibGit2.GitRepo, anns::Array{Base.LibGit2.GitAnnotated,1}, fastforward::Bool) at libgit2/merge.jl:96
  • merge!(repo::Base.LibGit2.GitRepo, anns::Array{Base.LibGit2.GitAnnotated,1}) at libgit2/merge.jl:80
  • mod(A::BitArray, B::Array{Bool,N<:Any}) at bitarray.jl:1132
  • mod(A::Array{Bool,N<:Any}, B::BitArray) at bitarray.jl:1133
  • next(a::Array, i) at array.jl:380
  • nextprod(a::Array{Int64,1}, x) at combinatorics.jl:169
  • permute!{Tv,Ti,Tp<:Integer,Tq<:Integer}(A::SparseMatrixCSC{Tv,Ti}, p::AbstractArray{Tp,1}, q::AbstractArray{Tq,1}, C::SparseMatrixCSC{Tv,Ti}, workcolptr::Array{Ti,1}) at sparse/sparsematrix.jl:944
  • permutedims!{T,N}(P::Array{T,N}, B::Union{Base.ReshapedArray{T,N,A<:DenseArray,MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N<:Any}}},DenseArray{T,N},SubArray{T,N,A<:Union{Base.ReshapedArray{T<:Any,N<:Any,A<:DenseArray,MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N<:Any}}},DenseArray},I<:Tuple{Vararg{Union{Base.AbstractCartesianIndex,Colon,Int64,Range{Int64}},N<:Any}},L<:Any}}, perm) at multidimensional.jl:983
  • prepend!(A::Array{Bool,1}, items::BitArray{1}) at bitarray.jl:693
  • prepend!{T}(a::Array{T,1}, items::AbstractArray{T<:Any,1}) at array.jl:499
  • process_exited(s::Array{Base.Process,1}) at process.jl:660
  • process_running(s::Array{Base.Process,1}) at process.jl:656
  • push!(a::Array{Any,1}, item::ANY<:Any) at array.jl:486
  • push!{T}(a::Array{T,1}, item) at array.jl:479
  • rand!{T<:Union{Bool,Int128,Int16,Int32,Int64,Int8,UInt128,UInt16,UInt32,UInt64,UInt8}}(rd::RandomDevice, A::Array{T,N<:Any}) at random.jl:53
  • rand!(r::MersenneTwister, A::Array{Float64,N<:Any}) at random.jl:348
  • rand!(r::MersenneTwister, A::Array{Float64,N<:Any}, n::Int64) at random.jl:348
  • rand!{I<:Base.Random.FloatInterval}(r::MersenneTwister, A::Array{Float64,N<:Any}, n::Int64, ::Type{I}) at random.jl:348
  • rand!{T<:Union{Float16,Float32}}(r::MersenneTwister, A::Array{T,N<:Any}, ::Type{Base.Random.Close1Open2}) at random.jl:375
  • rand!{T<:Union{Float16,Float32}}(r::MersenneTwister, A::Array{T,N<:Any}, ::Type{Base.Random.CloseOpen}) at random.jl:399
  • rand!{T<:Union{Float16,Float32}}(r::MersenneTwister, A::Array{T,N<:Any}) at random.jl:407
  • rand!(r::MersenneTwister, A::Array{UInt128,N<:Any}) at random.jl:411
  • rand!(r::MersenneTwister, A::Array{UInt128,N<:Any}, n::Int64) at random.jl:411
  • rand!{T<:Union{Int128,Int16,Int32,Int64,Int8,UInt16,UInt32,UInt64,UInt8}}(r::MersenneTwister, A::Array{T,N<:Any}) at random.jl:439
  • rand!(rng::MbedTLS.CtrDrbg, buf::Array) at /home/pedro/.julia/v0.5/MbedTLS/src/ctr_drbg.jl:40
  • read!(s::IO, a::Array{UInt8,N<:Any}) at io.jl:242
  • read!{T}(s::IO, a::Array{T,N<:Any}) at io.jl:247
  • readbytes!(ctx::MbedTLS.SSLContext, buf::Array{UInt8,1}) at /home/pedro/.julia/v0.5/MbedTLS/src/ssl.jl:208
  • readbytes!(ctx::MbedTLS.SSLContext, buf::Array{UInt8,1}, nbytes::UInt64) at /home/pedro/.julia/v0.5/MbedTLS/src/ssl.jl:210
  • readbytes!(s::IOStream, b::Array{UInt8,N<:Any}) at iostream.jl:232
  • readbytes!(s::IOStream, b::Array{UInt8,N<:Any}, nb) at iostream.jl:232
  • readbytes!(io::Base.AbstractIOBuffer, b::Array{UInt8,N<:Any}) at iobuffer.jl:326
  • readbytes!(io::Base.AbstractIOBuffer, b::Array{UInt8,N<:Any}, nb::Int64) at iobuffer.jl:328
  • readbytes!(io::Base.AbstractIOBuffer, b::Array{UInt8,N<:Any}, nb) at iobuffer.jl:326
  • readbytes!(s::Base.LibuvStream, a::Array{UInt8,1}) at stream.jl:714
  • readbytes!(s::Base.LibuvStream, a::Array{UInt8,1}, nb::Int64) at stream.jl:716
  • readbytes!(s::Base.LibuvStream, a::Array{UInt8,1}, nb) at stream.jl:714
  • readbytes!(ctx::MbedTLS.SSLContext, buf::Array{UInt8,1}, nbytes) at /home/pedro/.julia/v0.5/MbedTLS/src/ssl.jl:208
  • reinterpret{T,S}(::Type{T}, a::Array{S,1}) at array.jl:73
  • reinterpret{T,S}(::Type{T}, a::Array{S,N<:Any}) at array.jl:79
  • reinterpret{T,S,N}(::Type{T}, a::Array{S,N<:Any}, dims::Tuple{Vararg{Int64,N}}) at array.jl:86
  • reprmime(m::MIME{Symbol("text/vnd.graphviz")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME{Symbol("text/latex")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME{Symbol("text/calendar")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME{Symbol("text/n3")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME{Symbol("text/richtext")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME{Symbol("text/x-setext")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME{Symbol("text/sgml")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME{Symbol("text/tab-separated-values")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME{Symbol("text/x-vcalendar")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME{Symbol("text/x-vcard")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME{Symbol("text/cmd")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME{Symbol("text/css")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME{Symbol("text/csv")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME{Symbol("text/html")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME{Symbol("text/javascript")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME{Symbol("text/markdown")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME{Symbol("text/plain")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME{Symbol("text/vcard")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME{Symbol("text/xml")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME{Symbol("application/atom+xml")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME{Symbol("application/ecmascript")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME{Symbol("application/json")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME{Symbol("application/rdf+xml")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME{Symbol("application/rss+xml")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME{Symbol("application/xml-dtd")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME{Symbol("application/postscript")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME{Symbol("image/svg+xml")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME{Symbol("application/x-latex")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME{Symbol("application/xhtml+xml")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME{Symbol("application/javascript")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME{Symbol("application/xml")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME{Symbol("model/x3d+xml")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME{Symbol("model/x3d+vrml")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME{Symbol("model/vrml")}, x::Array{UInt8,1}) at multimedia.jl:57
  • reprmime(m::MIME, x::Array{UInt8,1}) at multimedia.jl:75
  • reshape{T,N}(a::Array{T,N}, dims::Tuple{Vararg{Int64,N}}) at array.jl:101
  • reshape{T,N}(a::Array{T,N<:Any}, dims::Tuple{Vararg{Int64,N}}) at array.jl:112
  • scale!{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64}}(X::Array{T,N<:Any}, s::T) at linalg/dense.jl:13
  • scale!{T<:Union{Complex{Float32},Complex{Float64}}}(X::Array{T,N<:Any}, s::Real) at linalg/dense.jl:27
  • scale!{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64}}(X::Array{T,N<:Any}, s::Number) at linalg/dense.jl:25
  • scale!{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64}}(s::T, X::Array{T,N<:Any}) at linalg/dense.jl:23
  • serialize(s::AbstractSerializer, a::Array) at serialize.jl:187
  • setindex!(A::Array{Any,N<:Any}, x::ANY<:Any, i::Int64) at essentials.jl:112
  • setindex!{T}(A::Array{T,N<:Any}, x, i1::Real) at array.jl:415
  • setindex!{T}(A::Array{T,N<:Any}, x, i1::Real, i2::Real, I::Real...) at array.jl:416
  • setindex!{T}(A::Array{T,N<:Any}, X::Array{T,N<:Any}, I::UnitRange{Int64}) at array.jl:444
  • setindex!(A::Array, X::AbstractArray, I::AbstractArray{Int64,1}) at array.jl:427
  • setindex!(A::Array, x, I::AbstractArray{Int64,1}) at array.jl:420
  • setindex!{T}(A::Array{T,N<:Any}, X::Array{T,N<:Any}, c::Colon) at array.jl:454
  • setindex!{T,N}(A::Array{T,N}, x::Number, ::Vararg{Colon,N}) at array.jl:463
  • setindex!(A::Array, x::Number, ::Colon) at array.jl:462
  • setindex!(B::BitArray, X::Union{Array,BitArray}, I0::Union{Colon,UnitRange{Int64}}) at multidimensional.jl:825
  • setindex!(B::BitArray, X::Union{Array,BitArray}, I0::Union{Colon,UnitRange{Int64}}, I::Union{Colon,Int64,UnitRange{Int64}}...) at multidimensional.jl:845
  • similar{T}(a::Array{T,1}) at array.jl:120
  • similar{T}(a::Array{T,2}) at array.jl:121
  • similar{T}(a::Array{T,1}, S::Type) at array.jl:122
  • similar{T}(a::Array{T,2}, S::Type) at array.jl:123
  • similar{N}(a::Array, T::Type, dims::Tuple{Vararg{Int64,N}}) at array.jl:125
  • similar{T}(a::Array{T,N<:Any}, m::Int64) at array.jl:124
  • similar{T,N}(a::Array{T,N<:Any}, dims::Tuple{Vararg{Int64,N}}) at array.jl:126
  • size(a::Array) at array.jl:20
  • size(a::Array, d) at array.jl:17
  • sizeof(a::Array) at array.jl:31
  • sort!{O<:Union{Base.Order.ForwardOrdering,Base.Order.ReverseOrdering{Base.Order.ForwardOrdering}},T<:Union{Float32,Float64}}(v::Array{Int64,1}, a::Base.Sort.Algorithm, o::Base.Order.Perm{O,Array{T,1}}) at sort.jl:625
  • srand(r::MersenneTwister, seed::Array{UInt32,1}) at random.jl:129
  • stacktrace(trace::Array{Ptr{Void},1}, c_funcs::Bool) at stacktraces.jl:147
  • stacktrace(trace::Array{Ptr{Void},1}) at stacktraces.jl:147
  • start(A::Array) at array.jl:379
  • startswith(a::Array{UInt8,1}, b::Array{UInt8,1}) at strings/util.jl:34
  • stringmime(m::MIME{Symbol("text/vnd.graphviz")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME{Symbol("text/latex")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME{Symbol("text/calendar")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME{Symbol("text/n3")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME{Symbol("text/richtext")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME{Symbol("text/x-setext")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME{Symbol("text/sgml")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME{Symbol("text/tab-separated-values")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME{Symbol("text/x-vcalendar")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME{Symbol("text/x-vcard")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME{Symbol("text/cmd")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME{Symbol("text/css")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME{Symbol("text/csv")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME{Symbol("text/html")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME{Symbol("text/javascript")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME{Symbol("text/markdown")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME{Symbol("text/plain")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME{Symbol("text/vcard")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME{Symbol("text/xml")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME{Symbol("application/atom+xml")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME{Symbol("application/ecmascript")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME{Symbol("application/json")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME{Symbol("application/rdf+xml")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME{Symbol("application/rss+xml")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME{Symbol("application/xml-dtd")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME{Symbol("application/postscript")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME{Symbol("image/svg+xml")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME{Symbol("application/x-latex")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME{Symbol("application/xhtml+xml")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME{Symbol("application/javascript")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME{Symbol("application/xml")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME{Symbol("model/x3d+xml")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME{Symbol("model/x3d+vrml")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME{Symbol("model/vrml")}, x::Array{UInt8,1}) at multimedia.jl:58
  • stringmime(m::MIME, x::Array{UInt8,1}) at multimedia.jl:77
  • success(procs::Array{Base.Process,1}) at process.jl:610
  • symperm{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, pinv::Array{Ti,1}) at deprecated.jl:765
  • trace{T}(A::Array{T,2}) at linalg/dense.jl:130
  • transcode{T<:Union{Int32,UInt32}}(::Type{T}, src::Array{UInt8,1}) at c.jl:273
  • transcode(::Type{UInt16}, src::Array{UInt8,1}) at c.jl:284
  • transcode(::Type{UInt8}, src::Array{UInt16,1}) at c.jl:334
  • transcode{T<:Union{Int32,UInt16,UInt32,UInt8}}(::Type{T}, src::Array{T,1}) at c.jl:271
  • transcode{S<:Union{Int32,UInt32}}(::Type{UInt8}, src::Array{S,1}) at c.jl:275
  • unsafe_copy!(dest::Array{Base.Pkg.Resolve.MaxSum.FieldValues.FieldValue,1}, doffs, src::Array{Base.Pkg.Resolve.MaxSum.FieldValues.FieldValue,1}, soffs, n) at pkg/resolve/fieldvalue.jl:72
  • unsafe_copy!{T}(dest::Array{T,N<:Any}, doffs, src::Array{T,N<:Any}, soffs, n) at array.jl:50
  • unsafe_copy!(dest::BitArray, doffs::Integer, src::Union{Array,BitArray}, soffs::Integer, n::Integer) at bitarray.jl:457
  • unshift!{T}(a::Array{T,1}, item) at array.jl:537
  • utf8(v::Array{UInt8,1}) at deprecated.jl:49
  • write(ctx::MbedTLS.MD, buf::Array{UInt8,1}) at /home/pedro/.julia/v0.5/MbedTLS/src/md.jl:138
  • write(s::IO, a::Array{UInt8,N<:Any}) at io.jl:175
  • write{T}(s::IO, a::Array{T,N<:Any}) at io.jl:179
  • subscribe(socket::ZMQ.Socket, filter::Union{AbstractString,Array}) at /home/pedro/.julia/v0.5/ZMQ/src/ZMQ.jl:225
  • unsubscribe(socket::ZMQ.Socket, filter::Union{AbstractString,Array}) at /home/pedro/.julia/v0.5/ZMQ/src/ZMQ.jl:225

Getting back to our example, we are ready to plot the trajectory of the particle. We will use an array to keep track of the positions of the particle over time.


In [25]:
using Plots  # Importing the library we will use.

h0 = 0.0  # Initial height (m).
v0 = 10.0  # Initial velocity (m/s).

endOfSimulation = 2.2  # We will stop the simulation at 2.2s.

posList = []  # The array to store the positions.

Δt = 0.1  # Time step

for i in 0:Δt:endOfSimulation
    
    push!(posList, h0)  # Storing the position in posList array.
    
    # Updating position and velocity:
    h0, v0 = h0 + v0*Δt, v0 - localGravity(g₀, rₑ, h0)*Δt
    
    if h0 < 0  
        # If the particle returns to the ground, we stop the simulation:
        endOfSimulation = i
        break
    end
    
end

# To plot the results we use the plot function from Plots:
plot(0:Δt:endOfSimulation, posList, xlabel="Time", ylabel="Height", w=3, label="")


INFO: Precompiling module Reexport.
INFO: Precompiling module FixedSizeArrays.
INFO: Precompiling module RecipesBase.
INFO: Precompiling module PlotUtils.
INFO: Precompiling module PlotThemes.
INFO: Precompiling module Showoff.
INFO: Precompiling module StatsBase.
INFO: Precompiling module NaNMath.
Out[25]:

The objects reaches a height of nearly $5.5\ m$ when launched with a velocity of $10\ m/s$.

Putting it all together

There is the concept of Escape Velocity, which is the velocity something needs to be launched from a planet's surface to escape its gravitational influence (and getting lost in space forever). So, what is the escape velocity of our model? We can try to estimate it, empirically. First, let us convert our simulation code into a function.


In [26]:
function simul(h0, v0, endOfSimulation)

    posList = []  # The array to store the positions.

    Δt = 0.1  # Time step

    for i = 0:Δt:endOfSimulation

        push!(posList, h0)  # Storing the position in posList array.

        # Updating position and velocity:
        h0, v0 = h0 + v0*Δt, v0 - localGravity(g₀, rₑ, h0)*Δt

        if h0 < 0  
            # If the particle returns to the ground we stop the simulation
            endOfSimulation = i
            break
        end
    end
    0:Δt:endOfSimulation, posList  # This is the return of our function
    
end


Out[26]:
simul (generic function with 1 method)

Now, we plot the trajectory of the particle for different initial velocities.


In [67]:
fig = plot()
for v0 in 11500:500:11500
    time, toPlot = simul(0, v0, 400000)
    plot!(time, toPlot, label=string(v0) * " m/s", w=3)
end
fig


Out[67]:

We can see from the plot that the escape velocity is bigger than 10500 $m/s$. The actual Earth escape velocity is estimated to be 11200 $m/s$. You can play with the values to see if you can get close to the actual value from our model.

In this Notebook, we explored some basic features of Julia, namely, variables and function declarations, while and for-loops and used the Plots library to visualize some output.